Interpolative Approach for Solving the Anderson Impurity Model 
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A rational representation for the self-energy is explored to interpolate the solution of the Anderson 
impurity model in general orbitally degenerate case. Several constrains such as the Friedel's sum 
rule, positions of the Hubbard bands as well as the value of quasiparticle residue are used to establish 
the equations for the coefficients of the interpolation. We employ two fast techniques, the slave- 
boson mean-field and the Hubbard I approximations to determine the functional dependence of 
the coefficients on doping, degeneracy and the strength of the interaction. The obtained spectral 
functions and self-energies are in good agreement with the results of numerically exact quantum 
Monte Carlo method. 
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I. INTRODUCTION 



There has been recent progress in understanding 
physics of strongly correlated electronic systems and their 
electronic structure near a localization-delocalization 
transition through the development of dynamical mean- 
field theory (DMFT)i. Merging this computation- 
ally tractable many-body technique with realistic local- 
density-approximation (LDA)S based electronic struc- 
ture calculations of strongly correlated solids is promis- 
ing due to its simplicity and correctness in both band 
and atomic limits. At present, much effort is being 
made in this direction including the developments of a 
LDA+DMFT method^, LDA++ approach^, combined 
GW and DMFT theory™, spectral density functional 
theorjiSi as well as applications to various systems such as 
Lar.^Sr^TiOgi, V 2 OA Fe and Ni&, Ce±a, PuiLi£, tran- 
sition metal oxides^, and many others. For a review, see 
Ref. El 

Such ab initio DMFT based self-consistent electronic 
structure algorithms should be able to explore all space 
of parameters where neither dopings nor even degeneracy 
itself is kept fixed as different states may appear close to 
the Fermi level during iterations towards self-consistency. 
This is crucial if one would like to calculate properties 
of realistic solid state system where bandwidth and the 
strength of the interaction is not known at the beginning. 
It is very different from the ideology of model Hamiltoni- 
ans where the input set of parameters defines the regime 
of correlations, and the corresponding many-body tech- 
niques may be applied afterwards. Realistic DMFT sim- 
ulations of material properties require fast scans of the 
entire parameter space to determine the interaction for 
a given doping, degeneracy and bandwidth via the solu- 
tion of the general multiorbital Anderson impurity model 
(AIM) 15 . Unfortunately, present approaches based on ei- 
ther non-crossing approximation (NCA) or iterative per- 



turbation theory (IPT) are unable to provide the solution 
to that problem due to a limited number of regimes where 
these methods can be applied 1 . The quantum Monte 
Carlo (QMC) techniquoi*i& is very accurate and can cope 
with multiorbital situation but not with multiplet inter- 
actions. Also its applicability so far has been limited ei- 
ther to a small number of orbitals or to unphysically large 
temperatures due to its computational cost. Recently 
some progress has been achieved using impurity solvers 
that improve upon the NCA approximationiLi&iii, but 
it has not been possible to retrieve Fermi liquid behav- 
ior at very low temperatures with these methods in the 
orbitally degenerate case. 

As universal impurity solvers have not yet being de- 
signed in the past we need to explore other possibilities, 
and this paper proposes interpolative approach for the 
self-energy in general multiorbital situation. We stress 
that this is not an attempt to develop an alternative 
method for solving the impurity problem, but follow- 
up of the ideology of LDA theory where approximations 
were designed by analytical fitsSS to the quantum Monte 
Carlo simulations for homogeneous electron gas2i. Nu- 
merically very expensive QMC calculations for the im- 
purity model display smooth self-energies at imaginary 
frequencies for a wide range of interactions and dopings, 
and it is therefore tempting to design such an interpola- 
tion. We also keep in mind that for many applications a 
high precision in reproducing the self-energies may not 
be required. One of such applications is, for example, 
the calculation of the total energylQ^lZ^ which, as well 
known from LDA based experience, may not be so sensi- 
tive to the details of the one-electron spectra. As a result, 
we expect that even crude evaluations of the self-energy 
shapes on imaginary axis may be sufficient for solving 
many realistic total energy problems, some of which have 
appeared alreadj*ii*«2»ii. Another point is a computa- 
tional efficiency and numerical stability. Bringing full 
self-consistent loops with respect to charge densities 11 
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and other spectral functions require many iterations to- 
wards the convergency which may not need too accurate 
frequency resolutions at every step. However, the proce- 
dure which solves the impurity model should smoothly 
connect various regions of the parameter space. This is 
a crucial point if one would like to have a numerically 
stable algorithm and our new interpolational approach 
ideally solves this problem. 

In the calculations of properties such as the low-energy 
spectroscopy and especially transport more delicate dis- 
tribution of spectral weight is taken place at low ener- 
gies, and the imaginary part of the analytically continued 
self-energy needs to be computed with a greater preci- 
sion. Here we expect that our obtained spectral functions 
should be used with care. Also, in a few well distinct 
regimes, such, e.g., as very near the Mott transition, the 
behavior maybe much more complicated and more dim- 
cult to interpolate. For the cases mentioned above exten- 
sions of the interpolative methods should be implemented 
and its beyond the scope of the present work. 

We can achieve a fast interpolative algorithm for the 
self-energy by utilizing a rational representation. The 
coefficients in this interpolation can be found by forcing 
the self-energy to obey several limits and constrains. For 
example, if infinite frequency (Hartree-Fock) limit, posi- 
tions of the Hubbard bands, low-frequency mass renor- 
malization z, mean number of particles n as well as the 
value of the self-energy at zero frequency S(0) are known 
from independent calculation, the set of interpolating co- 
efficients is well defined. In this work, we explore the 
slave-boson mean-field (SBMF) approac h 22 i 23 i 24 i 25 and 
the Hubbard I approximation 26 to determine the func- 
tional dependence of these coefficients upon doping, de- 
generacy and the strength of the interaction U. We ver- 
ify all trends produced by this interpolative procedure in 
the regimes of weak, intermediate and strong interactions 
and at various dopings conditions. These trends are com- 
pared with known analytical limits as well as against cal- 
culations using the quantum Monte Carlo method. Also, 
compared with QMC are self-energies and spectral func- 
tions on both imaginary and real axes for selective values 
of dopings. They indicate that the SBMF approach can 
predict such parameters of interpolation as n, E(0)and z 
with a good accuracy while the Hubbard I method fails in 
a number of regimes. However, the functional form of the 
atomic Green function which appears within Hubbard I 
can be used to determine positions of atomic satellites 
which helps to impose additional constraints on our pro- 
cedure. 

Giving an extraordinary computational speed of this 
approach we generally find a very satisfactory accuracy 
in comparisons with the numerically more accurate QMC 
calculations. If an increased accuracy is desired our 
method can be naturally extended by imposing more 
constraints and by implementing more refined impurity 
solvers other than the ones explored in this work. 

The paper is organized as follows. In Section II we dis- 
cuss rational interpolation for the self-energy and list the 



constraints. In Section III we discuss methods for solv- 
ing Anderson impurity model such as the slave-boson 
mean field and the Hubbard I approximations which can 
be used to find these constraints. A brief survey of the 
QMC method used to benchmark our algorithm is also 
given. We present numerical comparisons of SBMF and 
Hubbard I techniques against the QMC simulations for 
such quantities as quasiparticle residue and multiple oc- 
cupancies. In Section IV we report the results of the 
interpolative method and compare the obtained spectral 
functions with the QMC. In Section V we discuss pos- 
sible improvements of the algorithm. Section VI is the 
conclusion. 



II. INTERPOLATIVE APPROACH 

To be specific, we concentrate on the Anderson impu- 
rity Hamiltonian 

N N 
H = £fJ2faf a+ 9 U J2 U * n i3 + J2 EkaC £a Cka 



^[v-(k)/+ CkQ + y Q (k) c + Q / Q ], 



(1) 



describing the interaction of the impurity level e/ with 
bands of conduction electrons E\^ a via hybridization 
V a (k). U is the Coulomb repulsion between different 
orbitals in the /-band. Inspired by the success of the 
iterative perturbation theoryi, in order to solve the An- 
derson impurity model in general multiorbital case, we 
use a rational interpolative formula for the self-energy. 
This can be encoded into a form 
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The coefficients a m , 6 m , or, alternatively, the poles P, 
zeroes Zm and S(oo) in this equation are to be deter- 
mined. The form J5J can be also viewed as a continuous 
fraction expansion but continuous fraction representation 
will not be necessary for the description of the method. 

Our basic assumption is that only a well distinct set 
of poles in the rational representation is necessary to 
reproduce an overall frequency dependence of the self- 
energy. Extensive experience gained from solving Hub- 
bard and periodic Anderson model within DMFT at var- 
ious ratios of the on-site Coulomb interaction U to the 
bandwidth W shows the appearance of lower and up- 
per Hubbard bands as well as renormalized quasiparticle 
peak in the spectrum of one-electron excitationsi. 

It is clear that the Hubbard bands are damped atomic 
excitations and to the lowest order approximation appear 
as the positions of the poles of the atomic Green function. 
In the SU(N) symmetry case which is described by the 
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Hamiltonian QJ, these energies are numerated by the 
number of electrons occupying impurity level, i.e. E n — 
Cffi + \\Jn(n — 1), and the atomic Green function takes 
a simple functional form 



positions of poles and zeroes via recalculating probabil- 
ities X n which is equivalent to the famous Hubbard I 
approximation (discussed in more detail in the next Sec- 
tion) . 



JV-l 



G at (to) 



Cn 1 (Xn + X n+ i) 

lo + fj, — E n+ i + E n 



(3) 



where X n are the probabilities to find an atom in con- 
figuration with n electrons while combinatorial factor 
C^ _1 = n \^Zn-i)\ arr i ves due to equivalence of all 
states with n electrons in SU(N). 

We can represent the atomic Green function J3J using 
the rational representation J2Jl, i.e. 
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(G) (G) 

where P„ are all N atomic poles, while Z\ ' denote TV— 
1 zeroes with TV being the total degeneracy. The centers 
of the Hubbard bands are thus located at the atomic 
excitations P„ = E n — E n -\ — fi = e/ — fi — U(n — 1). 
Using standard expression for the atomic Green function 
G a t{u) = 1/[cj + /x — 6f — E a t (a;)], we arrive to a desired 
representation for the atomic self-energy 
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Using this functional form for finite A(u>) modifies the 



We now concentrate on the description of the quasi- 
particle peak which is present in metallic state of the 
system. For this an extra pole and zero have to be 
added in Eq To see this, let us consider the Hub- 
bard model for the SU(N) case where the local Green 
function can be written by the following Hilbert trans- 
form Gf(u>) = H[lu + /i — ej — E((x))]. If self-energy life- 
time effects arc ignored, the local spectral function be- 
comes Nf(uj) = D[uj + fj, — Cf — E(cj)] where D is the 
non-interacting density of states. The peaks of the spec- 
tral functions thus appear as zeroes in Eq. JSJ and in 
order to add the quasiparticle peak, one needs to add 
one extra zero (denoted hereafter as X) to the numer- 
ator in Eq. JSJ. To make the self-energy finite at 
lu — ► oo one has to also add one more pole (denoted 
hereafter as P{ ) which should appear in the denomi- 
nator or Eq. (JSJ- Furthermore, frequently the Hartree 
Fock value for the self-energy can be computed sepa- 
rately and it is desirable to have a parameter in the func- 
tional form J^Jl which will allow us to fix E(oo). An ob- 
vious candidate to be changed is that self-energy pole in 
(JSJ) which is closest to Li equal zero. Let us denote this 

and rewrite the denominator of J5j as 

N-2 



(S) 



parameter as P. 

{lu - P x (E) )(w - P 2 (S) ) nl w - Z « G) ] where the product is 

n=l 

now extended over all zeroes of the atomic Green func- 
tions except the one closest to zero and two extra poles 



P^> and P} > can control the width of the quasiparticle 
peak and E(oo). Thus, we arrive to the functional form 
for the self-energy 
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We are now ready to list all constrains of our interpola- 
te 

tive scheme. To fix the unknown coefficients X, P± ', 

P2 , Pn°^ , in Eq. @ and to write down the linear 
set of equations for the coefficients a m , b m in Eq. J2J. 
we use the following set of conditions. 

a) Hartree Fock value E(oo). In the limit to — ► oo the 
self-energy takes its Hartree-Fock form 

E(oo) = U(N - l)(n). (7) 

Mean level occupancy (n) is defined as a sum over all 



Matsubara frequencies for the Green's function, i.e. 

(n)=Tj2Gf(iu)e i » 0+ , (8) 

where 

Gf{u>) = — ^7 > v , y (9) 

defines the impurity Green function and A(cu) is the hy- 
bridization function. 



4 



b) Zero-frequency value S(0). The so called Friedel 
sum rule establishes the relation between the total den- 
sity and the real part of the self-energy at zero frequency 



1 1 

2 + n 



arctg 



;/ + jRE(zO+) +jRA(^0+) 
3A(iO+) 



2ni n ' dz 



(10) 



c) Quasiparticle mass renormalization value 
d^fS/dul^—Q. The slope of the self-energy at zero 
frequency controls the quasiparticle residue, z using the 
following relationship 



doj 



1 



(11) 



Formally, constraints (6) and (c) hold for zero tempera- 
ture only but we expect no significant deviations in many 
regions of parameters as long as we stay at low enough 
temperatures. The behavior may be more elaborated in 
the vicinity of Mott transition 27 . 
d) Positions of Hubbard bands. 

As we discussed, in order that the self-energy obeys the 
atomic limit and places the centers of the Hubbard bands 
at the positions of the atomic excitations, we demand 
that 



^ g) + m- £ / = s(^ g) ). 



(12) 



This condition fixes almost all self-energy zeroes ! in 

(s~t\ 

Eq.@ to the poles P n ■ However, it alone does not en- 
sure that the weight is correctly distributed among the 
Hubbard bands and that the very distant Hubbard bands 
disappear. For this to occur, distant poles of Green func- 
tion have to be canceled out by nearby zeroes of the 



(G) 

Green function. It is clear that each pole Pn far from 
the Fermi level has to be accompanied by a nearby zero 

(G) 

Z n in order the weight of the pole be small. Thus, the 
self-energy has poles at the position of Green's function 
zeroes which can be encoded into the constrain 



[E(ZP)]-i = 0. 



(13) 



We want to keep this property of the self-energy for finite 
A(w) and thus demand that self-energy diverges (when 
lifetime effects are kept, it only reaches a local maxi- 
mum) at the zeros of the functional form © of G a t{w). 
Note that the relationship 113|) holds (approximately) for 
frequency lo larger than the renormalized bandwidth zW. 

(G) 

Therefore the information about one Z n which lies close 
to lo = is omitted and replaced by the information 
about S(oo), S(0) and z as it is done by separating P± ' 
and P 2 in the denominator of Eq. ©. 

We can now write down a set of linear equations for 
all unknown coefficients in the expression (J2J). There is 
total 2M + 2 of parameters a m and b m ,m — 0, M, where 
we can always set 6q = 1. The conditions a),b),c) give 



ao = E(0), 

bo = 1, 

ai - E(0)6i = 1 - z 

o-m - 6m£(oo) = 0. 



(14) 
(15) 
(16) 
(17) 

According to condition d) we can use all N poles P n 

(G) (G) 

and N — 2 zeroes Z n ■ The zero Z n closest to u> = 
is dropped out. This brings additional 2N — 2 equations 
for the coefficients and makes M = N as the degree of 
the rational interpolation which are written below 



J 



N 



N 



m=0 
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N 
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(18) 
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Note that while M may be rather large, the actual 
number of poles contributing to the self-energy behavior 
is indeed very small. We can directly see this from Eq. 
© which uses all N poles P n G) fulfilling Eq. (T^J and 
uses N 



2 zeroes Z n 



directly related to N — 2 poles 
Pm' ■ Clearly, when the spectral weight of the atomic ex- 

(G) 

citation becomes small, the corresponding P„ ' becomes 



(G) 

close to Zn and the cancellation occurs. Therefore in 
realistic situations when only the upper and lower Hub- 
bard bands have significant spectral weight along with 
the quasiparticle peak, the actual degree of the polyno- 
mial expansion is either two or three. However, it is 
advantageous numerically and cheap computationally to 
keep all poles and zeroes in Eq. 10 because the formula 
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automatically distributes spectral weight over all existing 
Hubbard bands. 

In the limit when C/ — > the self-energy automatically 
translates to the non-interacting one. The atomic poles 
get close to each other but, most importantly, their spec- 
tral weight goes rapidly to zero as it gets accumulated 
within the quasiparticle band. 

In the Mott insulating regime, the conditions b) and 

(C) (G) 

c) drop out while all poles P„ and zeroes Z\ can be 
used to determine the interpolation. However, in this 

( G) 

regime it does not matter whether one of Z„ ' closest to 
lo = is dropped out or kept, since we can always replace 
this information by information about S(oo). Therefore 
the Mott transition can be studied without changing the 
constraints. 

We thus see that in the insulating case the self-energy 
correctly reproduces the well-known result of the Hub- 
bard I method where the Green function is computed 
after Eq. Jj|J with atomic self-energy. If the lifetime 

(G) (G) 

effects are computed, the parameters Pk and Z„ be- 
come complex and the Hubbard bands will acquire an 
additional bandwidth. This effect is evident from the 
simulations using various perturbative or QMC impurity 
solvers and can be naturally incorporated into the inter- 
polative formulas © or JfJJl. However, in practical imple- 
mentation below we will omit it for illustrative purposes. 

Let us now discuss the quality of interpolation from 
the perspective of the high-frequency behavior for the 
self-energy. The latter can be viewed^ as expansion in 
terms of the moments £( m ), i.e., 

00 £(m) 



Most important for us is to look at highest moments 
which are given by the Hartree Fock value, Eq. Q in- 
volving single occupancy matrix (n), as well as the first 
moment 

£« = [{N-l)(N~2)(nn) + (N-l)(n)-(N-l) 2 (n) 2 ]U 2 , 

(20) 

containing a double occupancy matrix (nn) . We see that 
the interpolation in part relies on the accuracy in com- 
puting multiple occupancies which are the functionals of 
both atomic excitations and the hybridization function. 
In this regard, using exact atomic Green function to find 
poles P^ and zeroes Z^ as part of the constrained pro- 
cedure may not be as accurate since it would assume the 
use of atomic multiple occupancies which do not carry 
information about A(w).On the other hand, we can also 
use only a functional form of the atomic Green function 
where the multiple occupancies are computed in a more 
rigorous manner. In the next Section we will show how 
this can be implemented using the SBMF multiple oc- 
cupancies which will be found to be in better agreement 
with the quantum Monte Carlo data. 



Note that the moments themselves can be used 

in establishing the constraints for interpolation coeffi- 
cients. This would involve independent evaluations of 
(n), (nn), (nun), etc. as well as various integrals involv- 
ing hybridization function A(w). However, we may run 
into ill-defined numerical problem since high-frequency 
information will be used to extract the low-frequency be- 
havior. Therefore, it is more advantageous numerically to 
use some poles and zeros of G at (u>) as given by condition 
d) above. 

We thus see that the interpolational scheme is defined 
completely once a prescription for obtaining parameters 
such as E(0),z, (n) as well as poles Pn , and zeroes 
Zn is given. For this purpose we will test two com- 
monly used methods: SBMF method due to Gutzwiller 22 
as described by Kotliar and Ruckenstein 23 and the well- 
known Hubbard I approximation^. We compare these 
results against more accurate but computationally de- 
manding quantum Monte Carlo calculations and estab- 
lish the procedure to extract all necessary parameters. 

Note that once the constraints such as z are computed 
from a given approximate method, some of the quanti- 
ties such as the total number of particles, (n), and the 
value of the self-energy at zero frequency, E(0), can be 
computed fully sclf-consistently. They can be compared 
with their non-self-consistent values. If the approximate 
scheme already provides a good approximation for (n) 
and satisfies the Friedel sum rule, the self-consistency 
can be avoided hence accelerating the calculation. Indeed 
we found that inclusion of the self-consistency improves 
the results only marginally except when we are in a close 
vicinity to the Mott transition but here we do not expect 
that our simple interpolative algorithm is very accurate. 

We now give the description of the approximate meth- 
ods for solving the impurity model and then present 
the comparisons of our interpolative procedure with the 
QMC calculations. 



III. METHODS FOR SOLVING IMPURITY 
MODEL 

A. Quantum Monte Carlo Method 

The quantum Monte Carlo method is a powerful and 
manifestly not perturbative approach in either interac- 
tion U or the bandwidth W . In the QMC method one in- 
troduces a Hubbard-Stratonovich field and averages over 
it using the Monte Carlo sampling. This is a controlled 
approximation using different expansion parameter, the 
size of the mesh for the imaginary time discretization. 
Unfortunately it is computationally very expensive as 
the number of time slices and the number of Hubbard- 
Stratonovich fields increases. Also the method works best 
at imaginary axis while analytical continuation is less ac- 
curate and has to be done with a great care. Extensive 
description of this method can be found in Ref. 01 We 
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will use this method to benchmark our calculations using 
approximate algorithms described later in this Section. 



B. Slave-Boson Mean Field Method 

A fast approach to solve a general impurity problem 
is the slave-boson method2iSi2^. At the mean field 
level, it gives the results similar to the famous Gutzwillcr 
approximation— However, it is improvable by perform- 
ing fluctuations around the saddle point. This approach 
is accurate as it has been shown recently to give the ex- 
act critical value of U in the large degeneracy limit at 
half-filling^. 

The main idea is to rewrite atomic states consisting 
of n electrons |'y 1 , ... ,7 n ), < n < N with help of a 
set of slave-bosons {^n 1 '"' 7 "}- In the following, we as- 
sume SU{N) symmetric case, i.e., equivalence between 
different states |7i)---,7 n ) for fixed n. Formulae corre- 
sponding to a more general crystal-field case are given in 
Appendix B. The creation operator of a physical electron 
is expressed via slave particles in the standard manner—. 
In order to recover the correct non-interacting limit at 
the mean-field level, the Bosc fields ip n can be considered 
as classical values found from minimizing a Lagrangian 
L{ip n } corresponding to the Hamiltonian Q). Two La- 
grange multipliers A and A should be introduced in this 
way, which correspond to the following two constrains: 

N 

EOl = l, (21) 

71=0 

N 

£ nC^l = G g (ico)e^ 0+ = n. (22) 

n—0 iul 

The numbers ij) n are similar to the probabilities X n 
discussed in connection to the formula for the atomic 
Green function ©. We thus see the physical meaning 
of the first constrain which is the sum of probabilities to 
find atom in any state is equal to one, and the second 
constrain gives the mean number of electrons coinciding 
with that found from G g (ui) = (ui — A — & 2 A(mj)) _:l . 
A combinatorial factor — n <^L n )\ arrives due to 
assumed equivalence of all states with n electrons. 

Minimization of L{i/j n } with respect to tp n leads us to 
the following set of equations to determine the quantities 

[E n +A~n\]i> n +nbTj2 A(iuj)G g (iuj)[LR^ n _ 1 +^ n bL 2 } + 

iuj 

+ (N - n)bT &(iuj)G g (iuj)[R 2 bi>n + LR^ n+1 ] = 0, 

(23) 

where b = RL^2 n _ 1 C n _iip n ip n _i, determines the 
mass renormalization, and the coefficients L — (1 — 




7T 

FIG. 1: Comparison between the slave boson mean field and 
the QMC calculation for (a) concentration versus chemical 
potential /i = n — e/ — (N — l)U/2, (b) dependence of the 
spectral weight Z on concentration, and (c) density-density 
correlation function, (nn) versus filling, n, in the two-band 
Hubbard model in SU(4) and U = 4 = 2W. 



En=i C»^l)-V\ R=(l- ElaC^r 1 ^)- 1 / 2 are 
normalization constants as in Refs. 23 24 E n = etn + 
Un(n — l)/2 is the total energy of the atom with n elec- 
trons in SU (N) approximation. 

Eq. H23|) . along with the constrains Q21[l. ill'L'l) consti- 
tute a set of non-linear equations which have to be solved 
iteratively. In practice, we consider Eq. I|23l) as an eigen- 
value problem with A being the eigenvalue and ip n being 
the eigenvectors of the matrix. The physical root corre- 
sponds to the lowest eigenvalue of A which gives a set of 
ip n determining the mass renormalization Z = b 2 . Since 
the matrix to be diagonalized depends non-linearly on il) n 
via the parameters L,R, and b and also on A, the solu- 
tion of the whole problem assumes the self-consistency: 
i) we build an initial approximation to ip n (for exam- 
ple the Hartree-Fock solution) and fix some A, ii) we 
solve eigenvalue problem and find new normalized tp n , 
iii) we mix new ip n with the old ones using the Broyden 
method~i and build new L, R, and b. Steps ii) and iii) 
are repeated until the self-consistency with respect to ip n 
is reached. During the iterations we also vary A to obey 
the constrains. The described procedure provides a sta- 
ble computational algorithm for solving AIM and gives 
us an access to the low-frequency Green's function and 
the self-energy of the problem via knowledge of the slope 
of S>E(iw) and the value RS(O) at zero frequency. 

The described slave-boson method gives the following 
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expression for the self-energy: 

E(w) = (l-b- 2 )uj-e f + Xb- 2 . (24) 

The impurity Green function Gj(u) in this limit is given 
by the expression 

G f (u) = b 2 G g (uj). (25) 

As an illustration, we now give the solution of Eq. ill'MI) 
for non-degenerate case (N = 2) and at the particle- 
hole symmetry point with e/ — fj, = — ^(N — 1). Con- 
sider a dynamical mean-field theory for the Hubbard 
model which reduces the problem to solving the impu- 
rity model subject to the self-consistency condition with 
respect to A(ui). Starting with the semicircular density 
of states (DOS), the self-consistency condition is given 
by Eq. 123(1 . We obtain the following simplifications: 
L = R = \/2, A = 0, ip = ip 2 ,b = k^ x i>2 and G g (u) = 
[ u - [W. ftfGgiu)]- 1 . The sum TJ2 iw A(ioj)G g (iuj) ap- 
peared in Eq. (|23() scales as Wa/2 with the constant a 
being the characteristic of a particular density of states 
and approximately equal to -0.2 in the semicircular DOS 
case. The self-consistent solution of Eq. (|23H is therefore 
possible and simply gives = 32 ^ Va + \- The Mott 
transition occurs when no sites with double occupancies 
can be found, i.e. when ip = ijj 2 = 0- A critical value 
of U c = 8W\a\. For a w -0.2, this gives U c w 1.6W 
and reproduces the result for U c2 — 1A9W known from 
the QMC calculation within a few percent accuracy. As 
degeneracy increases, critical U is shifted towards higher 
values 28 . From numerical calculations we obtained the 
following values of the critical interactions in the half- 
filled case U c w ZW for N = 6 (p-level), U c w 4.5VF for 
N = 10 (d-level), and U c 6W for N = 14 (/-level). 

Density-density correlation function (nn) for local 
states with n electrons is proportional to the number of 
pairs formed by n particles C 2 jC 2 . Since the probability 
for n electrons to be occupied is given by: P n = ip^C^ , 
the physical density-density correlator can be deduced 
from: (nn) = ^ n C 2 /C 2 P n . Similarly, the triple occu- 
pancy can be calculated from (nun) = C3 1 /C% P n . 

Let us now check the accuracy of this method by com- 
paring its results with the QMC data. We consider the 
two-band Hubbard model in SU(N — 4) orbitally degen- 
erate case. Hybridization A(cj) = V(\s.)/(u>— E^) sat- 
isfies the DMFT self-consistency condition of the Hub- 
bard model on a Bethe lattice 

A( W )=^j G( W ), (26) 

The Coulomb interaction is chosen to be U — 2W which 
is sufficiently large to open the Mott gap at integer 
fillings. All calculations are done for the temperature 
T = 1/32W. 

We first compare the average number of electrons 
vs. chemical potential determined from the slave bosons 
which is plotted in Fig. ^a). This quantity is sensitive 



to the low-frequency part of the Green function which 
should be described well by the present method. We see 
that it reproduces the QMC data with a very high accu- 
racy and only differs by 20 per cent very near the jump 
of fj at n = 1. 

The quasiparticle residue z versus filling n is plotted in 
Fig. mb). The slave-boson method gives the Fermi liq- 
uid and provides estimates for the quasiparticle residue 
with the overall discrepancy of the order of 30%. In 
fact, we have performed several additional calculations 
for other degeneracies (N = 2 and 6) and for various pa- 
rameters regimes. The trend to overestimate mass renor- 
malization can be seen in many cases. It disappears only 
when U approaches zero. We need to point out, how- 
ever, that i) the extractions of zero frequency self-energy 
slopes from the high-temperature QMC is by itself nu- 
merically not well grounded procedure, as information for 
the self-energy is known at the Matsubara points only, 
which is then extrapolated to w = 0, ii) other methods for 
solving impurity model, such as NCA or IPT display sim- 
ilar discrepancies and iii) recent findings 28 suggest that 
at least at half-filling quasiparticle residues deduced from 
slave bosons becomes exact when N — > 00. Most im- 
portantly for our interpolative method is that the entire 
functional dependence of z vs. filling, interaction and 
degeneracy is correctly captured. Its overall accuracy is 
acceptable as it is evident from our comparisons of the 
spectral functions presented in the next Section and well 
within the main goal of our work to provide fast scans of 
the entire parameter space necessary for simulating real 
materials. This is important as, for example, for general 
/-electron material, the QMC method is prohibitively 
time consuming, but we expect from the SBMF method 
the results for mass renormalization not worse than 50% 
for such delicate regime as the vicinity of the Mott tran- 
sition. 

Fig. ^c) shows the density-density correlation func- 
tion (nn) as a function of average occupation n. The 
discrepancy is most pronounced for fillings h < 1 [see 
the inset of Fig. ^c)] where the absolute values of (nn) 
are rather small. Although our slave-boson technique 
captures only the quasiparticle peak, it gives the correla- 
tion function in reasonable agreement with the QMC for 
dopings not too close to the Mott transition. 

C. Hubbard I Approximation 

Now we turn to the Hubbard I approximation 2 ^ which 
is closely related to the moments expansion method 30 . 
Consider many-body atomic states |$1 ) which in 
SU(N) are all degenerate with index n numerating these 
states for a given number of electrons n. The impurity 
Green function is defined as the average 

G f (r) = -(T T f a (r)f+(0)). (27) 

and becomes diagonal with all equal elements in 
SU(N)At is convenient to introduce the Hubbard opera- 
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tors The impurity Green function (|27() is given by 

Z$ = \*£WP\ (28) 

and represent the one-electron creation and destruction 
operators as follows 

/« = EE(*«°l/a|$i" +1) )^ +1 I (29) 
n kk' 

fi = EEf^ 1 '!/^!" 1 )^ 1 " (so) 

n kk' 

i 

where the matrix G nm (r) is defined as 

G nm (r) = - (^Val^ 15 )^^ 1 ^^^))^ 1 ^^!*^)^ (32) 



G f (r) =^G„(r), (31) 

n m 



Establishing the equations for G nm (r) can be per- 
formed using the method of equations of motion for 
the X operators. Performing their decoupling due to 
HubbardS^i, carrying out the Fourier transformation 
and analytical continuation to the real frequency axis, 
and summing over n and m after (|31|l we arrive to the 
net result 

Gj 1 (u) = G- t 1 (co)-A(u), (33) 

The Gati^j) can be viewed in the matrix form l|31|l 
with the following definition of a diagonal atomic Green 
function 

r at (,,\_z Cn~ 1 ( X n+ X n+1 ) 
LO + fl — £j n +l + M/ n 

with E n — efn + Un(n — l)/2 being the total ener- 
gies of the atom with n electrons in SU(N). The co- 
efficients X n are the probabilities to find atom with n 
electrons and were already discussed in connection to the 
formula J3J for the atomic Green function. They are sim- 
ilar to the coefficients ip^ introduced within the SBMF 
method but now found from different set of equations. 
These numbers are normalized to unity, J2n=o X n = 
J2n=o ^n^ 1 (^n + X n +i) — 1, and are expressed via di- 
agonal elements of G nm {iuj) as follows: 




FIG. 2: Comparison between the Hubbard I and the QMC 
calculation for (a) concentration versus chemical potential 
jx = jLt — ej — (N — l)U/2, (b) dependence of the spectral 
weight Z on concentration, and (c) density-density correla- 
tion function, (n a n a i) versus filling, n, in the two-band Hub- 
bard model in SU(A) and for U = 2W = 4. 



X n = -T^G n „(^)e-- 0+ /C _1 - (35) 

Their determination in principle assumes solving a 
non-linear set of equations while determining Gf(u>). 
The mean number of electrons can be measured as 
follows: n — ^2^=0 n ^n X n or as follows n = 
TNY,iu> G f {iuj)e^ a+ . The numbers X n can be also used 
to find the averages (nn), (nnn) in the way similar to 
what has been done in the SBMF approach. 



If we neglect by the hybridization A(w) in Eq. 133fl . 
the probabilities X n become simply statistical weights: 

e -E n /T 

X n = S? ■ (36) 

We thus see that in principle there are several different 
ways to determine the coefficients X n , either via self- 
consistent determination (|35|l . or using statistical for- 
mula l|36|) . or taking them from SBMF equation 123|l . 
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i.e. setting X n = i\) 2 n while still utilizing the functional 
dependence provided by the Hubbard I method. To de- 
termine the best procedure let us first consider limits 
of large and small C/'s. When A(w) = 0, Gf(ui) is 
reduced to ^ nm G^ m {u)), i.e. the Hubbard I method 
reproduces the atomic limit. Setting [7 = gives 
Gf(w) = [uj + n — tf — A((jj)] _1 , which is the correct 
band limit. Unfortunately, at half-filling this limit has 
a pathology connected to the instability towards Mott 
transition at any interaction strength U. To see this, we 
consider a dynamical mean-field theory for the Hubbard 
model. Using semicircular density of states, we obtain 
G f (u) = [1-(^) 2 G/ {u)G at {u)]- x G at {u) and conclude 
that for any small U the system opens a pathological 
gap in the spectrum. Clearly, using Hubbard I only, the 
behavior of the Green function at to — > cannot be re- 
produced and the quality of the numbers X n is at ques- 
tion. This already emphasizes the importance of using 
the slave-boson treatment at small frequencies. 

Ultimately, making the comparisons with the QMC 
calculations is the best option in picking the most accu- 
rate procedure to compute the probabilities X n . To check 
the accuracy against the QMC we again consider the two- 
band Hubbard model in SU (4) symmetry as above. The 
chemical potential, mass renormalization and double oc- 
cupancy are plotted versus filling in Fig. [21 All quan- 
tities here were computed with statistical weights after 
Eq. (|36|l but we found a similar agreement while utilizing 
the self-consistent determination of X n after Eq. (|35fl . 
We first see that the Hubbard I approximation does not 
give satisfactory agreement with the QMC data for n(fx) 
because it misses the correct behavior at low frequencies. 

The comparisons for z(n) plotted in Fig. ^b) sur- 
prisingly show a relatively good behavior. However, the 
pathology of this approximation at half-filling would pre- 
dict z — for any U, which is a serious warning not to 
use it for extracting the quasiparticle weight. Fig. |2fc) 
shows (nn) as a function of average occupation n. As 
this quantity is directly related to the high-frequency 
expansion one may expect a better accuracy here. How- 
ever, comparing Fig. |5[c) and Fig. [Jc), ^ is clear that 
the slave boson method gives more accurate double oc- 
cupancy. This is due to the fact that the density matrix 
obtained by the slave boson method is of higher quality 
than the one obtained from the Hubbard I approxima- 
tion. 

The results of these comparisons suggest that the prob- 
abilities ip^ provided by the slave-boson method is a bet- 
ter way in determining the coefficients X n in the metal- 
lic region of parameters. Therefore it is preferable to 
use these numbers while establishing the equations for 
the unknown coefficients in the interpolational form J2Jl. 
However, the functional form l|34l) of the Hubbard I ap- 
proximation with X„ = can still be used as it pro- 

(G) (G) 

vides the positions of the poles Pn and zeroes Z) x of 
the atomic Green function necessary for the condition 
d) in the previous Section. This also ensures accurate 
high-frequency behavior of the interpolated self-energy 



since its moments expressed via multiple occupancies are 
directly related to ip n . 

Interestingly, while more sophisticated QMC approach 
captures both the quasiparticle peak and the Hubbard 
bands this is not the case for the slave-boson mean field 
method. To obtain the Hubbard bands in this method 
fluctuations need to be computed, which would be very 
tedious in a general multiorbital situation. However the 
slave-boson method delivers many parameters in a good 
agreement with the QMC results, and, hence, it can be 
used to give a functional dependence of the coefficients 
of the rational approximation. 



IV. RESULTS OF THE INTERPOLATIVE 
SCHEME 

By now the procedure to determine the coefficients is 
well established. We use the SBMF method to deter- 
mine n, £(oo), S(0), z as well as poles and zeroes of the 
atomic Green function provided by the SBMF probabili- 
ties 4>n a n d by the bare atomic energy levels E n (we omit 
the lifetime effects for simplicity) . This generates a set of 
linear equations for coefficients , 6^ in the rational in- 
terpolation formula J3J. In the present Section we show 
the trends our interpolative algorithm gives for the spec- 
tral functions in various regions of parameters as well as 
provide detailed comparisons for some values of doping 
for both imaginary and real axis spectral functions. The 
two-band Hubbard model with semicircular density of 
states and DMFT self-consistency condition after l|2t)|) is 
utilized in SU(N = 4) symmetry in all cases using the 
bandwidth W — 2 and temperature T = 1/16. 



A. Trends 

Fig. 13 shows the behavior of the density of states 
N(u>) = —^sGf(oj)/ir for U = W as a function of 
the chemical potential jx computed with respect to the 
particle-hole symmetry point [N — l)U/2 and as a func- 
tion of frequency uj. The semicircular quasiparticle band 
is seen at the central part of the figure. Its band- 
width is only weakly renormalized by the interactions 
in this regime. It is half-filled for fi = (i.e. when 
/i = ef — (N — l)U/2) and gets fully emptied when chem- 
ical potential is shifted to negative values. Several weak 
satellites can be also seen on this figure which are due to 
atomic poles. Their spectral weight is extremely small 
in this case and any sizable lifetime effect (which is not 
included while plotting this figure) will smear these satel- 
lites out almost completely. While approaching fully 
emptied (or fully filled situation) the spectral weight 
of the Hubbard bands disappears completely and only 
unrenormalized quasiparticle band remains. It is clear 
that even without shifting the atomic poles to the com- 
plex axis, the numerical procedure of generating the self- 
energy is absolutely stable. 
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FIG. 3: Calculated density of states using the interpolative 
method as a function of chemical potential fi = y — ef — 
(N — l)U/2 and frequency for the two-band Hubbard model 
in SU(4) and at U = W = 2. 
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FIG. 4: Calculated density of states using the QMC method 
as a function of chemical potential fj, = y, — e/ — (TV — 1)?7 /2 
and frequency for the two-band Hubbard model in SU(4) and 
at U = W = 2. 
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FIG. 6: Calculated density of states using the QMC method 
as a function of chemical potential p, = y, — e/ — (AT — l)U/2 
and frequency for the two-band Hubbard model in SU(4) and 
at t/ = 2W = 4. 
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This trend can be directly compared with the simula- 
tions using a more accurate QMC impurity solver. We 
present this in Fig. 01 for U = W, which shows calcu- 
lated density of states in the same region of parameters. 
Remarkably that again we can distinguish the renormal- 
ized quasiparticle band and very weak Hubbard satellites. 
The Hubbard bands appear to be much more diffuse in 
this figure mainly due to the lifetime effects and partially 
due to maximum entropy method using for analytical 
continuation from imaginary to real axis. Otherwise the 
entire picture looks very much like the one on Fig. 
generated with much less computational effort. 

Fig. El gives the same behavior of the density of states 
for the strongly correlated regime U = 2W. In this case 
the situation at integer filling is totally different as the 
system undergoes metal-insulator transition. This is 
seen around the dopings levels with fl between and 
-1 and between -3 and -5 where the wight of the quasi- 
particle band collapses while lower and upper Hubbard 
bands acquire all spectral weight. In the remaining region 
of parameters both strongly renormalized quasiparticle 
band and Hubbard satellites remain. Again, once full 
filling or full emptying is approached the quasiparticle 
bands restores its original bandwidth while the Hubbard 
bands disappear. The QMC result for the same region 
of parameters is given in Fig. Again we can distin- 
guish the renormalized quasiparticle band and Hubbard 
satellites as well as the areas of Mott insulator and of 
strongly correlated metal. The Hubbard bands appear to 
be more sharp in this figure which signals on approaching 
the atomic limit. 



B. Comparison for Spectral Functions at Imaginary 
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FIG. 7: Comparison between real and imaginary parts of the 
Green function obtained from the interpoaltive method and 
the quantum Monte Carlo calculation for the two-band Hub- 
bard model at a set of fillings n = 0.5,0.8,1.2,1.5,1.8 and 
U = 2W = 4. 



We now turn to the comparison of the Green functions 
and the self-energies obtained using the formulae 10, and 
(ffl) respectively against the predictions of the quantum 
Monte Carlo method. We will report our comparisons 
for the two-band Hubbard model and sets of dopings 
corresponding to n — 0.5, 0.8, 1.2, 1.5, 1.8 using the value 
of U = 2W = 4. Other tests for different degeneracies, 
doping levels and the interactions have been performed 
which display similar accuracy. 

Fig. [7| shows the comparison between the real and 
imaginary parts of the Green function obtained by the 
interpolative method with the results of the QMC cal- 
culations. As one can see almost complete agreement 
has been obtained for a wide regime of dopings. The 
agreement gets less accurate once the half-filling is ap- 
proached, but still very good giving an extraordinary 
computational speed of the given method compared to 
QMC. 

Fig. |51 shows similar comparison between the real and 
imaginary parts of the self-energies obtained by the inter- 
polative and the QMC method. We can see that the self- 
energies exhibit some noise which is intrinsic to stochas- 



tic QMC procedure. The values of the self-energies near 
iuj = and ilu = oo are correctly captured with some 
residual discrepancies are attributed to slightly differ- 
ent chemical potentials used to reproduce given filling 
within every method. The results at the imaginary axis 
show slightly underestimated slopes of the self-energies 
within the interpolative algorithm which is attributed to 
the underestimated values of z obtained from the SBMF 
calculation. Ultimately improving these numbers by in- 
clusions of fluctuations beyond mean field will further 
improve the comparisons. However, even at the present 
stage of the accuracy all functional dependence brought 
by the SBMF method quantitatively captures the behav- 
ior of the self-energy seen from time consuming QMC 
simulation. 



C. Comparison for Spectral Functions at Real Axis 

We also made detailed comparisons between calculated 
densities of states obtained at the real axis using the in- 
terpolative method and the QMC algorithm. The QMC 
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FIG. 8: Comparison between real and imaginary parts of the 
self-energies obtained from the interpoaltive method and the 
quantum Monte Carlo calculation for the two-band Hubbard 
model at a set of fillings n = 0.5,0.8, 1.2, 1.5, 1.8 and U = 
2W = 4. 
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FIG. 9: Comparison between the one-electron densities of 
states obtained from interpolative formula and the quantum 
Monte Carlo calculation for the two-band Hubbard model at 
fillings n = 0.5, 0.8, 1.2, 1.5, 1.8. and U = 2W = 4. 



densities of states require an analytical continuation from 
the imaginary to real axis and were generated using the 
maximum entropy method. By itself this procedure in- 
troduces some errors within the QMC especially at higher 
frequencies. In Fig. EI we show our results for the fillings 
corresponding to n = 0.5, 0.8, 1.2, 1.5, 1.8 using the value 
of U = 2W = 4. One can see the appearance of the quasi- 
particle band and two Hubbard bands distanced by the 
value of U. It can be seen that the interpolative method 
remarkably reproduces the trend in shifting the Hubbard 
bands upon changing the doping. It automatically holds 
the distance between them to the value of U while this 
is not always true in the quantum Monte Carlo method. 
Despite this result, the overall agreement between both 
methods is very satisfactory. 

V. DISCUSSION 

Here we would like to discuss possible ways to further 
improve the accuracy of the method. The inaccuracies 
are mainly seen in three different quantities: i) the width 
of the Hubbard bands, ii) the mass renormalization z{fx) 



which is borrowed from the SBMF method, and iii) the 
number of electrons n(/i) extracted from the interpolated 
impurity Green function @ . The inaccuracy in the width 
of the Hubbard band is mainly connected to neglecting 
the lifetime effect. Provided it is computed, this will 
shift the positions of atomic poles onto the complex plane 
which is in principle trivial to account for within our in- 
terpolative algorithm. To improve the accuracy of z(jjl) 
one can, for example, work out a modified slave boson 
scheme which will account for the fluctuations around 
mean field solution. The inaccuracy in n(ji) is small in 
many regions of parameters and typically amounts to 5— 
10 per cent. We can try to improve this agreement by re- 
quirement that n(/z) obtained via interpolation matches 
with nsBMF{n) obtained by the SBMF method. The 
latter agrees very well with the QMC for a wide region 
of parameters as it is evident from Fig. Ufa). In reality, 
our analysis shows that in many cases the discrepancy in 
n(/x) is connected with the overestimation of zsbmf{^)- 
Therefore, points ii) and iii) mentioned above are inter- 
related. 

The requirement that = nsBMF^) can be en- 

forced by adjusting the width of the quasiparticle band, 
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and in many regions of parameters this is controlled 
by z(p). However, there are situations when the Hub- 
bard band appears in the vicinity of u> = 0, and chang- 
ing z(n) does not affect the bandwidth. To gain a con- 
trol in those cases it is better to replace the constraints 
S(0), dY, / dLu\uj = o by constraints of fixing the self-energies 
at two frequencies, say 2(0) and S(iwo) where loq is the 
frequency of the order of renormalized bandwidth. We 
have found that this scheme brings mass renormalizations 
which are about 30% smaller than the SBMF ones, and 
the agreement with the QMC is significantly improved. 
Thus, inaccuracies ii) and iii) can be avoided with this 
very cheap trick. However, we also would like to point 
out that the condition n{jj) — nsBMF(n) is essentially 
non-linear as the solution may not exist for all regions of 
parameters. It is, for example, evident that in such points 
where n(/i) is given by a symmetry (as, e.g., particle-hole 
symmetry point n = 2 in the case considered above) the 
mass renormalization does not affect the number of elec- 
trons. 

As the philosophy of our approach is to get the best 
possible fit we are also open to implementing any kinds 
of ad hoc renormalizations constants. One of such possi- 
bility could be the use of a quasiparticle residue 30 per 
cent smaller than zsbmf(/J>)- As z(fi) should go to unity 
when U = 0, the correction can, for example, be encoded 
into the formula z(fj,) — zsbmf((J-)[0-7 + 0.3zsbmf(h)}- 

We finally would like to remark that the scheme defined 
by a set of linear equations for the coefficients (|14[l 119|) 
is absolutely robust as solutions exist for all regimes of 
parameters such as strength of the interaction, doping 
and degeneracy. In general, bringing any information on 
the self energy £(w x ) at some frequency point u) x or its 
derivative d£/<iw| w=Wx would generate a linear relation- 
ship between the interpolation coefficients, thus keeping 



robustness of the method. On the other hand, fixing 
such relationships as numbers of electrons brings non- 
linearity to the problem which could lead to multiplicity 
or non-existence of the solutions. It is also clear that by 
narrowing the regime of parameters, the accuracy of the 
interpolative algorithm can be systematically increased. 



VI. CONCLUSION 

To summarize, this paper shows the possibility to inter- 
polate the self-energies for a whole range of dopings, de- 
generacies and the interactions using a computationally 
efficient algorithm. The parameters of the interpolation 
are obtained from a set of constraints in the slave boson 
mean field method combined with the functional form 
of the atomic Green function. The interpolative method 
reproduces all trends in remarkable agreement with such 
sophisticated and numerically accurate impurity solver as 
the QMC method. We also obtain a very good quantita- 
tive agreement in a whole range of parameters for such 
quantities as mean level occupancies, spectral functions 
and self-energies. Some residual discrepancies remain 
which can be corrected provided better algorithms deliv- 
ering the constraints will be utilized. Nevertheless, given 
the superior speed of the present approach, we have ob- 
tained a truly exceptional accuracy times efficiency of the 
proposed procedure. 
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VII. APPENDIX 

In the crystal field case we assume that AT— fold degenerate impurity level e/ is split by a crystal field onto G sublevels 
e/i, —ef a , ...e/G- We assume that for each sublevel there is still some partial degeneracy d a so that J2 a =i da = N. 
In limiting case of SU(N) degeneracy, G = l,c?i = N, and in non-degenerate case, G = N,di = d a = da = 1. We 
need to discuss how a number of electrons n can be accommodated over different sublevels e f a . Introducing numbers 
of electrons on each sublevel, n a , we obtain X* =i n a — n. Note the restrictions: < n < N, and < n a < d a . In 
SU(N) case, G = 1, n\ = n, and in non-degenerate case, G = N , n a is either or 1. Total energy for the shell with 
n electrons depends on particular configuration {n a } 

G 1 

E ni ...n G = ^ e /" n « + ^UCS a n a )[CS a n a ) - 1]. (37) 

a=l 

Many-body wave function is also characterized by a set of numbers {n a }, i.e. \nx...nG)- Energy E ni _.. na remains 
degenerate, which can be calculated as product of how many combinations exists to accommodate electrons in each 
sublevel, i.e. x ...C^° x ...C^° . Let us further introduce probabilities ip ni ,,, nG to find a shell in a given state with 
energy E ni ,„ na . Sum of all probabilities should be equal to 1, i.e. 



di d a dc 



(38) 
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There are two Green functions in Gutzwiller method: impurity Green function Gf{iu>) and quasiparticle Green 
function G g {ui) — b Gf(td)b , where matrix coefficients b represent generalized mass renormalizations parameters. 
All matrices are assumed to be diagonal and have diagonal elements enumerated as follows: Gi(ui), ...G a {u>), ...Gg{oj). 
Each element in the Green function is represented as follows 



1 

w - X a - b%A a (u)' 



G ga (u) = — ^r^, (39) 



G fa (w) = b 2 a G ga {u>). (40) 
and determines a mean number of electrons in each sublevel 

n a = d a TY,G ga {iuy M+ . (41) 

iuj 

The total mean number of electrons is thus: n — }] h a . Hybridization function A(iw) is the matrix which 
is assumed to be diagonal, and it has diagonal elements enumerated as follows: Ai(w), ...A a (uj), ...Ag(w). Mass 
renormalizations Z a = b 2 a are determined in each sublevel. 



Here: 



-enerj 


ry are 




= w + 


IJ,-ef a 


- A tt (w) 


di 






<E 


... 5: .. 


• E c t 


»1=0 


n a = l 


n G =0 


di 




da 


E ■ 


- E ■■■ 


E 


ni=0 


n a = l 


n G =0 


di 


do-1 


<i G 


E- 




E 


ni=0 


n Q =0 


n G =0 



1 . Aq, 

52"^ e / Q- 52"' 



a=l 
G 



(42) 



(43) 

. -1/2 



(44) 



\ -1/2 

R a = I 1- 2J • E ^-^nr 1 - C 'nS<...n a ... nG • (45) 

\ ni— n a — n G — / 

The generalization of the non-linear equations (|23|1 has the form 
= [E ni ... nG + A- (Y,%\ a n a )] Vni...no + 

G 

^n Q [TS iw A Q (iw)G 9Q ,(icj)]6 a [-R Q! iaV'„ 1 ...n c «-i...„ G + b aL 2 a ip ni .. ,„ q ... „ g ] + 



^(d Q - n a )[TS M A ct (iw)G gQ (iw)]6 Q V'n 1 ...n a +l...r iG + ^ai?aV'„ 1 ...n ...n G ] ■ ( 4 6) 

Q = l 
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